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Abstract 

The molecular dynamics and Monte Carlo studies of the thermal stability of 
C36, and Ceo fullerenes and the methane molecule are reported. It has been 
shown that the heat transfer between the atomic cluster and the external heat 
reservoir can either promote or prohibit the formation of defects in this cluster. 
The wide temperature and pressure ranges have been determined where the defect 
formation rate is essentially non Arrhenius. It has been shown that the lifetime 
of light clusters in molecules depends more strongly on the contact with the heat 
reservoir. A statistical model that is based on the kinetic equation and allows 
for an analytical solution has been developed to explain the results. Within this 
model, the generalized Arrhenius formula has been derived to predict the lifetime 
of the clusters in an arbitrary thermal contact with the environment. 
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The mean lifetime r of an atomic cluster to its decay or transition to another 
^ ! isomer is often described by the Arrhenius law 

g : r(T) = £e E «/ kT , (1) 

where k —is the Boltzmann constant, E a — is the decay activation (isomerization) 
'k! ' energy, T — is the temperature of the environment in contact with the cluster, 
and the preexponential factor £ has the dimension of time and is the inverse 
of the frequency factor of the corresponding process. In some approaches, this 
preexponential factor is considered constant (see, e.g., [T]), but in the general 
case it depends not only on the atomic cluster structure and the characteristic 
heat transfer time Tth between the cluster and environment, but also on the 
temperature, i.e., £ = £(7^, T) (see, e.g., [21 E])- 

If there is no thermal contact between the cluster and the environment {jth 
+00), the temperature dependence of r has another fundamentally nonArrhenius 
form |H 

^=(f-tr- « e >j- ( 2) 

[ +00, if E < E a 

Here, E = CT m , E = Eki n +E pot — is the sum of the kinetic and potential energies 



of the cluster (the latter energy is measured from the energy of the equilibrium 
state of the cluster), where C — is the heat capacity of the cluster and T m — is 
the microcanonical temperature [6], [7], relating to the average kinetic energy of 
the cluster atoms (Eki n ) as (Ekin) = \kT m X, where X = (3N at — 6) — is the 
number of degrees of freedom of the cluster consisting of N at atoms taking into 
account the conservation of three components of the total momentum and the 
angular momentum. 

As the heat transfer intensity between the cluster and environment increases, 
lifetime dependence (ED, is transformed into Eq. (QQ). The key objective of this 
work is to study the character of the (ED to (DQ) transformation, as the heat transfer 
time Tth decreases to a finite value, and to determine the relation between the 
preexponential factors £ and £ of the decay of the thermalized and heatinsulated 
clusters. 

In a computer simulation, the quantity r is determined as the average of 
the lifetimes of several clusters decayed at the same temperature (the energy 
for the heatinsulated clusters) but at different initial distributions of the atomic 
velocities in these clusters. If the number of clusters Mq is sufficiently large, then 
the number of clusters M(i), retained by the time t, decreases exponentially as 
M(t) = Moexp(—t/r). At a sufficiently large number of statistical tests, the 
probability density n finding the cluster in a certain state can be described in the 
framework of the kinetic equation [H]. 

For the cluster of N at atoms, the n value is generally a function of 6N at 
coordinates of the phase space and the time [8J. To solve the main problem 

of this work, it suffices to assume that n depends only on the time t and the 

+00 

total configuration energy E, i.e., M(t) = J n(E,t)dE. Deriving the kinetic 



equation, it is covenient to consider that the clusters are introduced into the 
system continuously at a constant rate I rather than simultaneously at the time 
t = 0. In this case, the initial energy of the cluster is Eq = CT (or CT m for the 
heatinsulated cluster). It is also assumed that the cluster lifetime is much longer 
than the characteristic period of its natural oscillations and the characteristic time 
of the internal heat transfer between its normal modes. The kinetic equation of 
the particle balance is 

/ +00 
= IS(E-E ) I_n(£,t)+— l-n(E,t) + S(E,T) [ n(E\t)dE> 

Ot Tq{E) T th \ J 

. (3) 

where Eq = CT, S(E, T) — is the equilibrium distribution density of an individual 
thermalized cluster in the energy E at a temperature T. To avoid the consider 
ation of the heat transfer processes from the environ ment to the individual atoms 
of the cluster, we assumed that like in [D], the entire cluster is simultaneously 
thermalized; i.e., the nonequilibrium distribution of its states with the density 



n(E, t) is transformed in the characteristic time interval Tth to the equilibrium 
one with the density S(E, T). 

The dependence of the average cluster lifetime r on the characteristic thermaliza- 
tion time Tth an d other parameters of interest can be obtained from the stationary 
solution no(E) of Eq. 



-oo 



= I5(E-E ) L^ no (E,t) + —l-n (E) + S(E,T) [ n (E')dE' . (4) 

ME) Tth \ J J 

For the steadystate process, the average lifetime r, the total number of the clusters 
M(t) \t=+oo in the system, and the rate I of their addition to the system are related 
as 

+00 

J n (E)dE = rI. (5) 


Substituting Eq. (J5j) into Eq. (HI), we obtain 



r g(£,r)r + r t „a(E-E ) 
ME) = 7 1 + r th /r (E) ' (6) 

Integrating expression @ with respect to E and taking into account Eq. (j5J) and 

+00 

the relation J S(E,T)dE = 1, we obtain the desired dependence of r on the 


parameters of the problem in the form 

+00 



The preexponential factor £ is independent of the activation energy but can 
depend on other parameters, such as the temperature [3], atomic mass, characteris- 
tic thermalization time r^, etc. In the frequent collision limit, — > +00), 
the cluster decay can be considered as the diffusion of the atoms with a certain 
diffusion coefficient D from the metastable equilibrium point through the energy 
barrier separating the cluster from its decay products [21 HD]- In this case, the 
lifetime depends on D. Meanwhile, Eq. (J7j) in the limit l/r t h — > +00 has the form 

+00 

l/r= J S(E,T)(l/r (E))dE. (8) 


In this expression, the only parameter that can be a function of is the preexpo- 
nential factor £ in Eq. ([2]). For further analysis, it suffices to expand this parameter 
into the Taylor series up to the linear term in 1/D. Taking into account that the 
diffusion coefficient for Brownian motion is proportional to the mean free path 



time and temperature [8], and the mean free path time is proportional to the 
thermalization time r^, we have D ~ t^T, and the preexponential factor £ can 
be represented as 

e«G> + 6/(7toD. (9) 

We studied the transition of the heatinsulated decay mode to the thermalized 
one on an example of the four clusters of C6o, C36, C20 fullerene and CH4 methane 
molecule. For the Ceo and C36 fullerenes, we studied the isomerization process 
resulting from the Stone- Wales transformation [H] . The decay process was conside- 
red for the C20 fullerene and CH4 molecule. The C20 fullerene is the lightest of the 
existing fullerenes; for this reason, methane is added to the series of the clusters 
under study in order to follow the effect of the cluster mass on the character of 
the transition from law (J2D to law (PQ). 

The fullerenes were simulated with the orthogonal tight-binding potential [12], 
successfully used in the simulation of both small clusters and macroscopic carbon 
structures (earlier, we applied it to investigate various characteristics of the 
and C20 fullerenes and oter carbon systems JI31 (HJ d5l dSl fTTj). To determine the 
lifetime of the CH4 molecule, we used the nonorthogonal tightbinding method |18| , 
which requires much more computer resources than the orthogonal tight-binding 
method, but was successfully used to determine the activation energy of the 
hydrocarbon structures |T§] . 

The parameters E a: £0 an d £1, in Eq. (J7j), were determined by the numerical 
simulation of the decay or isomerization of the corresponding clusters. The activa- 
tion energy E a of the processes of interest was determined by the static simulation 
method. Initially, the stationary points of the interaction potential in the (3N at — 
6)-dimensional space of the generalized coordinates of the cluster atoms were 
found. These stationary points correspond to the metastable and saddle configura- 
tions whose energies are used to determine the height of the energy barrier 
separating the initial configuration from the decay products (or the nearest isomer 
if the isomerization process is considered). The process activation energy was 
assumed as equal to this barrier height. The details of the static simulation method 
are presented in [20]. In the general case, the activation energy can differ from the 
barrier height, but these two quantities for fullerenes coincide with each other |22| . 
We also observed a similar pattern in methane. 

We determined the cluster lifetime at different temperatures and heattransfer 
intensities by the molecular dynamics (MD) method and/or by the combined 
molecular dynamics and Monte Carlo (MDMC) method. The MD method was 
used for the temperatures at which the cluster lifetime did not exceed several 
nanoseconds. In this case, the Newton equation describing the motion of cluster 
atoms was numerically solved by the velocity Verlet method [22] The heat transfer 
with the environment was simulated by random collisions of the cluster atoms with 
the buffer gas atoms as in [D]. In this work, we were interested in a wide range of the 
cluster thermalization times r t h from 10~ 16 s to 1 s. The fullerenes are synthesized 



at the buffer gas (helium or argon) pressure, which provides r t h ~ 10~ 8 s [23] , and 
the fast thermalization with the characteristic time r t h = 10~ 14 — 10~ 16 s can be 
obtained at the cluster thermal contact with the crystalline substrate. The MD 
method allows one to determine the lifetime (or the isomerization time) of the 
cluster at an arbitrary (including zero) intensity of the heat transfer between the 
cluster and the environment, but it requires large computer resources and is thus 
applicable only for sufficiently high temperatures at which the lifetime is short 
(~ Ins). The MDMC method was applied to numerically determine the cluster 
lifetime at low temperatures at which the lifetime is ~ Is [21]; this method 
is applicable at a sufficiently high heat transfer intensity when the Tth value is 
comparable with the minimum period of the natural oscillations of a fullerene that 
is 24 and 20 fs for C20 and Ceo respectively, in the tightbinding approximation. 
In their combination, the MD and MDMC methods allow one to investigate the 
entire ranges of temperatures and heattransfer intensities. 

To compare the numerically simulated average cluster lifetimes to Eq. (J7J), the 
cluster heat capacities C and the distribution densities S(E, T) were determined 
in the harmonic approximation 

S(E ' T) = X'.ikT)^ * 6 '^' ° = kX - (10) 



Approximation ( TTOjl is justified by the data reported in [21], where it was shown 
that the energy distribution density in the C20, fullerene calculated in the tight- 
binding model is almost identical to Eq. ( TTQT ). 

The parameters E a , £0 and £1 for the C20, C36 and Ceo fullerenes and the 
methane CH4 molecule are presented in the table. The activation energy E a was 
determined by the static simulation method to a relative accuracy e ~ 10~ 4 . The 
£0 value was determined both by the comparison of Eq. (J2J) and the lifetime for 
the heatinsulated clusters obtained by the MD method and by the comparison of 
Eq. (J7j) and the lifetime for the thermalized clusters obtained by the MDMC and 
MD methods. The £1 value was determined by comparing relation (J7j) and the 
lifetime for the thermalized clusters obtained by the MDMC and MD methods. 
The number of numerical calculations N var , determining the statistic error of £0 
and £1, is also presented in the table. 

The tabulated data imply that the £, value determined by the MD method in the 
heatinsulated clusters coincides to a statistical error with the £0, value determined 
by the MDMC method for the clusters in contact with the heat reservoir. This 
indicates both the validity of the MDMC approach and the adequacy of analytical 
dependence (J7j). The accuracy of the determination of £0 by the MDMC method is 
higher than that for the MD method, because the MDMC method covers a wider 
temperature range. 

We note that the dependences of the lifetimes of the methane molecule and 
fullerenes on the heattransfer intensity are different. Relation (J7j) implies that the 
lifetime of the heatinsulated cluster at sufficiently low temperatures is infinite, 



Table 1: Activation energies E a , preexponential factor £ and parameters £ an d £1 of the decay (isomerization) 
preexponential factor for the C20, C36. and Ceo fullerenes and methane molecule 





E a ,eV 
Method 


(1/T t fc=0) 

Method (N var ) 


to, s 
(1/rtfc ^0) 
Method (N var ) 


a, s 2 K 
(l/r t fc ^0) 
Method (N var ) 


C60 


6.48 
SM 


jq— 17.31±0.91 

MD(40) 


jq— 17.17±0.26 

MDMC(16) 


jq— 3b.l3±1.17 

MDMC(16) 


C36 


5.74 
SM 


^0-15.72=^0. 82 
MD(40) 


jQ-15.58iO.19 
MDMC(24) 


jQ-35.33i0.87 
MDMC(24) 


C20 


5.00 
SM 


jQ-17.60iO.83 

MD(40) 


jq-17.46±0.14 

MDMC(24) 


jq-35.2I±0.79 

MDMC(24) 


CH 4 


3.55 
SM 


MD(100) 


jQ-16.34iO.02 

MD(4000) 


jq-34.6S>±0.43 

MD(4000) 



r(rth,T)\i/ Tth= Q = 00 at T < T c , T c = E a /kX (since the cluster is heat insulated, 
the microcanonical temperature is assumed). For the methane molecule and C20, 
C 36 , and C 60 fullerene, T c = 4576 K, 1074 K, 653 K and 432 K, respectively. Therefo- 
re, at the typical temperatures of petrochemical reactions T ~ 800 K, the depen- 
dence of the lifetime of the methane molecule on its heat transfer intensity is 
always nonmonotonic. Figure 1 presents this dependence at different temperatures. 
The figure implies that the methane molecule is least stable at the characteristic 
heat transfer time with the environment r t h ~ 10~ 13 — 10~ 14 s at any temperature. 

The heattransfer intensity dependence of the lifetime of comparatively heavy 
fullerenes has a different shape. Figure 2 presents the r(T t h, T) curves for the 
PF20 fullerene. The comparison of Figs. 1 and 2 shows that the decrease in the 
lifetime for the PF20 fullerene with increasing heatexchange intensity is not so 
large, the minimum is more shallow and is similar to a plateau; i.e., the fullerene 
lifetime is independent of Tth in a wide range of this parameter. For the C36 
and Ceo fullerenes, this trend gets more pronounced, because the clusters with a 
larger number of atoms are themselves heat reservoirs, and the defect formation 
processes in them depend much weaker on the interaction with the environment. 

To conclude, we note that the dependence of the cluster lifetime on the heat 
transfer conditions presented in this work is more universal than the previous 
approximations. However, its applicability field is also limited by the assumption 
of the constant heat capacity of the cluster and its quasiharmonic behavior to 
the decay time. It would be interesting to extend this approach to the case of a 
molecule with pronounced anharmonic oscillations. 

We are grateful to L.A. Openov for useful discussions and critical comments. 
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Fig.l. Log-log plot of the lifetime r of the methane molecule versus the charac- 
teristic heattransfer time rth with the thermostat at different temperatures. The 
curves are the theoretical dependence (7) for temperatures (from top to bottom) 
10 000, 7000, 4000, 3500, 3000, and 2000 K. The circles, squares and rhombuses 
are the numerical calculation data for 4000, 3500, and 3000 K, respectively. Each 
symbol is the result of averaging over 100 independent calculations. 




Fig. 2. Log-log plot of the lifetime r of the C20 fullerene versus the characteristic 
heattransfer time r t h with the thermostat at different temperatures. The curves 
are the theoretical dependence (7) at temperatures (from top to bottom) 4000, 
2824, 2182, 1777, 1500, and 1300 K. The circles are the numerical calculation 
data. Each circle is the result of averaging over 16 independent calculations. 



